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ABSTRACT 

Aims. We investigate the interplay between the ionization radiation from massive stars and the turbulence inside the surrounding 
molecular gas thanks to 3D numerical simulations. 

Methods. We used the 3D hydrodynamical code HERACLES to model an initial turbulent medium that is ionized and heated by an 
ionizing source. Three different simulations are performed with different mean Mach numbers (1,2 and 4). A non-equilibrium model 
for the ionization and the associated thermal processes was used. This revealed to be crucial when turbulent ram pressure is of the 
same order as the ionized-gas pressure. 

Results. The density structures initiated by the turbulence cause local curvatures of the dense shell formed by the ionization com- 
pression. When the curvature of the shell is sufficient, the shell collapse on itself to form a pillar while a smaller curvature leads 
to the formation of dense clumps that are accelerated with the shell and therefore remain in the shell during the simulation. When 
the turbulent ram pressure of the cold gas is sufficient to balance the ionized-gas pressure, some dense-gas bubbles have enough 
kinetic energy to penetrate inside the ionized medium, forming cometary globules. This suggests a direct relation in the observations 
between the presence of globules and the relative importance of the turbulence compared to the ionized-gas pressure. The probability 
density functions present a double peak structure when the turbulence is low relative to the ionized-gas pressure. This could be used 
in observations as an indication of the turbulence inside molecular clouds. 
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of clum ps studied the formation of globules and lMackey & LinJ 
(|201Q|) proposed the shadowing effect to explain the formation 
of pillars. 

In a previous paper (iTremblin et alJ (l2012h . paper I here- 
after), we presented a new approach emphasizing the importance 
of the curvature of the dense shell at the edge of the HII region 
to explain the formation of pillars and dense clumps in the shell. 
When the shell is sufficiently curved, a pillar will form by the 
collapse of the shell on itself, while a smaller curvature will trig- 
ger lateral flows forming dense clumps and dips inside the shell. 
However this scenario and the previous ones presented above are 
rather idealized set-ups and have to be validated in more realistic 
situations, e.g. taking i nto account a turbulent molecular cloud. 

R ecent studies (see | Mellema et al.ll20()6l: iGritschneder et al.l 

l201Ql:lArthur et al.l[20T ll) started to concentrate on the interplay 
between ionization from the massive stars and the turbulence in- 
side the molecular cloud. They found that pillars, dense clumps 
and globules emerge naturally in their models, however the de- 
tailed processes forming these structures are difficult to identify 
with the turbulence. 

In the present study, we concentrate on the interplay between 
turbulence and ionization in the light of the detailed processes 
presented in paper I. We first present the numerical method used 
and some tests of the out-of-equilibrium ionization module, fol- 
lowed by three ionized turbulent simulations at different mean 
Mach numbers (1,2 and 4) and finally the comparison with the 
previous idealized situations presented in paper I and a new ap- 
proach to explain the formation of globules. Finally we show 



1. Introduction 

Although massive stars have a great impact on their envi- 
ronment, the importance of their radiative feedback on star- 
formation ratesi still a matter of discussion. For example, 
iDale et al.l (l2QQ5h found a slight enhancement of the star- 
formation activity when the radiative feedback is taken into 
account whereas they recently presented simulations on a 
scale of a giant mole cular cloud with almost no impact (see 
iDale & Bonnellll20TTI) . A closer look at the small scales helps 
addressing this question by identifying the detailed mechanisms 
that form the dense structures eventually leading to star forma- 
tion. These structures consist mainly of pillars of gas pointing 
toward the ionizing source, globules detached from the parent 
molecular cloud, and dense clumps at the interface between the 
HII region and the cloud. Optical, infrared (IR), and far-IR ob- 
servations, using the Hubble space telescope and the Spitzer 
and Herschel satellites impressively revealed th ese features. 



mainly in high-mass s tar-forming regions (e.g. iHesteretal 



19961: lOerin et al.ll2009l: IPeharveng et al.ll2QlQl: [Schneider et al 



2Q10HZavagno et al Several models have been proposed 



to explain their formation. The collect and collapse scenario de- 
scribed by[Elmegreen & Lada (19 ^ or sha dowing instabilities 
in the ionization front (e.g. IWilliamslll999l) concentrate on the 
formation of dense clump s at the interface. iBertoldil (Il989l) and 
iLefloch & Lazarefj (1 19941) with the radiation driven implosion 
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that probability density functions of the gas can be used in ob- 
servations to determine the importance of turbulence relative to 
the ionized-gas pressure. 

2. Numerical methods 



We used the HERACLES cod43 (iGonzalez et al.ll2007h with the 
ionization n iodule described in paper I and the co oling module 
described in I Audit & HennebelQ (l20ld) (based on lWolfire et alJ 
(l2003h ) to model the molecular cloud and the ionization coming 
from the OB cluster. In addition to the physics of paper I, the 
turbulence in t he cloud is modeled usin g the turbulence mod - 
ule described in lSchmidt etaP (l2006h andi lSchmidt et al.1 (l2009l) . 
Gra vity is an importan t ingredient for ultra-compact HII region 
(e.g. lPeters et al.ll20TQl) . however we are interest ed in large scale 
HII re gions (e.g. the Rosette Molecular cloud, ISchneider et al.l 
(120101) . Self-gravity and the gravity from the ionizing source 
were considered in paper I without noticeable change on the for- 
mation of the structures. Therefore we do not consider gravity in 
the present study. 

First, we describe the model we used for ionization. It is an 
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Fig. 1. Cut of the density field and velocity field in the radiation 
driven implosion of a clump inside an ionized medium. The gas 
behind the clump in the white box is out of equilibrium and re- 
combines slowly. The gas at the front of the clump is escaping 
by the rocket eff'ect and is cooled down by the strong expansion. 

out-of-equilibrium model solving the equations 

d(Fy)/dx = -rinil - X)cryFy -h ^o(-^)^* 
d(XnH)/dt = I-R = nH(l- X)(TyFy - pX^n]j 

de/dt = Ixey-Rx hT/(y - 1) (1) 

in which X is the ionization fraction (advected with the hy- 
drodynamic), Fy is the photon flux, the source term, hh 
the total hydrogen density, CTy the cross-section for ionization, 
/ the ionization rate, R the recombination rate, is given by 
2 X 10"^^r"^-^^cm^/s, e is the internal energy of the gas, Cy 
the energy transferred from the ionizing photons to the elec- 
trons and T the temperature of thermal equilibrium between all 
the species. The first equation is the ray tracing of the ionizing 
photons coming from the source. The second one is the ioniza- 
tion/recombination photo-chemistry evolution and the third one 
is the thermodynamic evolution associated. ID and 3D tests of 
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Fig. 2. ID and 3D simulations of an isolated clump inside a 
completely ionized medium exposed to an ionization flux of 10^ 
photons per second per squared centimeter. Note that in the ID 
simulation the clump is in fact an infinite slab. Top: ionization 
profile through the middle of the box at t = 500 ky. Bottom 
Temperature profile at T = 500ky. The ionization source is at 
X = parsec. In the ID simulation the clump is located at x = 
1.5 parsec while in the 3D simulations it is at 2.75 parsec. The 
two 3D simulations show that it is not needed to limit too much 
the variation of the ionization fraction. The gas behind the clump 
is in the shadow and is therefore being recombined. Its recombi- 
nation time is long and is not at equilibrium. These simulations 
also show the importance of the 3D eff'ects on the temperature 
structure. 



the out-of-equilibrium ionization and recombination processes 
are performed in Sect. [3] 

The turbulent forcing is modeled by an Ornstein-Uhlenbeck 
process in Fourier space. For each mode of the force field, ran- 
dom increments with a Gaussian distribution are added after 
each time step. The wave-numbers of the forcing modes are in 
the interval [0, 2ko], where ko is given by the typical length scale 
of the forcing, L = In/ko with L = 4 pc, the size of the box in 
our case. By projecting the modes in transversal and longitudinal 
directions with regard to the wave vectors, divergence-free and 
rotation-free components are produced in physical space. The 
resulting force field is statistically homogeneous, isotropic, and 
stationary. The RMS magnitude is of the order V^/L, where V is 
the characteristic velocity of the flow. We will perform three sim- 
ulations with V resulting in stationary transonic and supersonic 
velocity: Mach 1, 2 and 4. The turbulence used is solenoidal with 
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a ratio of compressive forcing power to the total forcing power 
of 10 %. 

Three turbulent scenarii are investigated in Sect. IH a tran- 
sonic turbulent molecular cloud (mean Mach number around 1) 
and two supersonic cases (mean Mach number of 2 and 4). The 
box is a cube of 4 parsec at a resolution of 0.01 pc (400^), the 
mean density is 500 cm"^ and the temperature is initially at 25 
Kelvin. During the turbulent evolution, the boundary conditions 
are periodic. The simulations are ran until a statistically steady 
state is obtained, i.e. a constant mean Mach number in the box. 
Then we turn on the ionization with a flux of 10^ photons per 
second per squared centimeter, which is the typical flux of an 
04 star at 30 parsecs. The ionization is coming from the top of 
the box and this face is changed to reflexive boundary condi- 
tion for the hydrodynamic while the opposite side is set to free 
flow. The turbulent stirring is kept during the ionization phase. 
Some runs were done without the turbulence maintained with- 
out noticeable change, the large-scale modes do not have time to 
influence the small scales during the ionization phase. In Sect.lH 
we investigate the interplay between turbulence and ionization 
for the three difl'erent simulations. 

3. Tests of the out-of-equilibrium ionization and 
recombination processes 

When turbulence is included in an ionization simulation, the 
main difl'erence is the possibility that some ionized gas gets into 
the shadow of the dense unionized gas because of the mixing. 
This is typically an out-of-equilibrium state for the ionized gas 
that will begin to recombine. To study this state, we investigate 
simplified set-ups in which hot ionized gas is put in the shadow 
of dense cold gas and therefore will be out of equilibrium. 
We simulate the impact of ionization on a clump (t)ciump=0.5 
parsec) inside an homogeneous ionized medium (see Fig.[l]). We 
also perform an "equivalent" ID simulation where the clump 
will be an infinite slab, to identify the eff'ect of the geometry. 
The gas in the shadow of the clump/slab will start to recombine 
since it is out of equilibrium. The typical recombination time 
for the ionized gas is given by I /fin that is of the order of 
10^-10^ years for ionized gas at a density of 1-10/cm^. It is of 
the order of the gas dynamical time: Ddump/Cio = 10^ years, 
in which C/^ is the sound speed of the ionized gas. Therefore 
an out-of-equilibrium model is needed to treat this state of the 
gas. We studied first the profile of the ionization fraction and 
the eff'ect of the limited- variation time- stepping and then the 
structure of the temperature in these tests. 

Our method solves the two first equations in ([T]) implicitly 
while the thermal balance is solved explicitly and is sub-cycled 
to limit the variation of T. Usually the variation of X is limited 
at 10 % in the implicit step. However, it is often not needed 
to limit that much the variations. The ionization-fraction and 
the temperature profiles are plotted in Fig. [2l after 500 ky 
of evolution and for limited variations of 10 and 500 % and 
there is almost no diff'erence between the two. At the front, the 
equilibrium for ionization is reached even with long time- step 
thanks to the implicit method, the ionization fraction is nearly at 
one. At the back, the recombination time of a hot ionized plasma 
is quite long. For example in our case the recombination time is 
of the order of 100 ky which also explains why the ionization 
fraction is still around 0.2 after 500 ky. Therefore it is possible 
to allow long time-steps. 

We now investigate the temperature structure. The 3D 
profiles in Fig. [2] are taken through the middle of the box (the 
white-dashed line in Fig [T]), therefore through the middle of 



the clump, which is located between 2.5 and 3 parsecs. In the 
ID case, the clump is a slab located after 500 ky between 1 
and 2 parsecs. The diff'erence of position between the two is 
easy to understand: in the ID case, the ionized material cannot 
be evacuated on the side, the column density in front of the 
ionization front is higher. There are only a few ionization 
processes occurring at the ionization front, the dynamics is 
dominated by the expansion of the gas evaporating from the 
slab. In the 3D case, the ionized gas is escaping on the side of 
the clumps, therefore ionization can penetrate further inside the 
clump and deposit more energy. As a consequence, the clump is 
pushed further away by rocket motion and the energy deposit at 
the surface of the clump leads to a high expansion of the ionized 
gas and to a shock surface that can be seen in Fig.[T] The ionized 
gas is cooled down to 3000-6000 Kelvin by the expansion at the 
front. 

The temperature structure is also quite diff'erent behind 
the clump (white box in Fig. [T]) because of the shadowed 
recombination. In the ID simulation, the thermal equilibrium is 
reached at the front, in the HII region.Behind the slab, the gas 
is not at equilibrium and recombines, the ionization fraction is 
decreasing from 1 to 0.2. The pressure drops by nearly a factor 
of 2 when the ionization fraction decreases from 1 to ^ 0.2. 
Moreover, most of the energy emitted during the recombination 
process is radiated away by the cooling, so that the temperature 
of the gas does not vary much, increasing only slightly, from ^ 
8000 K to ^ 9000 K. The temperature structure is quite diff'erent 
in 3D as shown in Fig. [J] The cooling by recombination behind 
the clump should decrease the pressure of the gas as in the ID 
case. However in 3D, it leads to lateral gas flows from the hot 
ionized gas exposed to the radiation to the hot recombined gas 
in the shadow. Therefore the recombination behind the clump is 
done more or less at constant pressure leading to an increase in 
temperature behind the clump, which is surprising at first sight. 

This simple test shows that the equilibrium assumption for 
ionization is invalid when hot ionized gas gets in the shadow 
of dense unionized gas (Ionization fraction of 20 % instead of 
0). The equilibrium assumption is usually justified by the fact 
that the recombination time is much shorter than the dynamical 
time. This is usually true in the HII region, however this is 
not the case at the interface between the HII region and the 
cold gas when turbulence is included. Indeed especially at high 
level of turbulence, hot ionized gas will be mixed with cold 
dense gas and therefore will be in the shadow of the cold gas, 
exactly in the situation we studied in this simplified test. The 
ionization equilibrium is not reached for a non-negligible part 
of the ionized gas that can get in the shadow of dense unionized 
gas. This is especially the case when the level of turbulence 
is sufficient to balance the ionized-gas pressure and to mix the 
unionized and ionized phases as we will see in the next section. 



4. Turbulent simulation 

4.1. Transonic turbulence 

The first line of Fig. [3] presents three snapshots of the column 
density for the transonic simulation, when ionization is turned 
on, 500 ky and 1 My after ionization has started. The hot ion- 
ized gas expands in the cold turbulent medium and triggers the 
propagation of a shock ahead of the ionization front. This is 
the collect p art of the cl assical collect and collapse scenario 
( Elmesreen & Ladall 19771) and is observed in n iany regions (see 
iZavagno et al.1 120101: [Thompson et al.l 120111: iDeharveng et al.l 



3 



Tremblin et aL: Globules and Pillars Simulations 




4.0 I 

3.5 

u 3.0! 

B 2.5 



^ I 2.01 

S 1.5 1 
Fi 1.0 B 



^ 0.5 



3 ^'^0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 



^ 0.5 
a 

•S ^-%.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 
^ [parsec] 



4.0 
3.5 

^ 3.0 
i-H 

2 5 
g ^ 2.0 



Time: My 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 





0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 





0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 
[parsec] [parsec] 

Time: 0.5 My Time: 1 My 



22.5 

22.4 

22.3 " 
B 

22.2 ^ 
22.1 

o 

22.0 Q 

1 ^ 

21.9 g 

I 21.8 S 
be 

I 21.7 a 

21.6 
' 21.5 



Fig. 3. Column density snapshots of the three turbulent simulations (row 1 : Mach 1, row 2: Mach 2, row 3: Mach 4) et three different 
time (column 1: before ionization starts, column 2: 500 ky after ionization is turned on, column 3: 1 My). 



I2OIOI) . The 500-ky snapshot clearly shows a three-phase medium 
with the hot ionized gas at the top, the shocked dense region and 
the unperturbed cold cloud. The shocked region is apparently 
0.5-1 parsec wide at 500 ky and gets wider in time, up to 1.5 
parsec at 1 My. This width is a lot larger than one would expect 
from the collect process in a homogeneous medium. The typical 
width of the dense shell is 0.1 parsec in this case. However Fig. 
|4] clearly shows that the real thickness of the shell is indeed of 
order 0.1 parsec but contrary to the collect process in a homoge- 
neous medium, the shell is not flat but highly curved. Therefore 
the column density plots show a wide shocked region just by 
projection efl'ects of this highly turbulent shocked region. 

The properties of these phases can be studied using the dis- 
tribution of the mass fraction in the density-pressure plane or 
in the density-Mach number plane. Fig. [5] and [6] show the 2D 
plots of the mass fraction before and after ionization is turned 
on. Before ionization is turned on, the distribution of the gas 



is at thermal equilibrium and the average Mach number is of 
order 1 with a dispersion of 0.48. After 1 My, the three-phase 
medium is formed and the three phases are well separated in the 
Mach number density plot. The hot ionized gas can be seen on 
the pressure density plot, on the isothermal equilibrium curve 
for the ionized gas, and on the Mach number density plot in the 
sub- sonic area. The unperturbed medium remains at the same 
location in the graph with respect to Fig. [5] while the shocked 
material is pushed at high Mach number and high densities. The 
distribution of the mass fraction is quite similar to the one found 
in paper I, the only difl'erence is the initial distribution of the gas, 
which was previously homogeneous at 500 cm"^ with a density 
or interface modulation. In both cases, the shocked region con- 
tains points with densities higher than the one corresponding to 
the collect and collapse scenario (red dashed line in Fig.[5]and[6]), 
i.e. the density after the collect phase in a homogeneous medium. 
These points represent 3 % of the shell in mass fraction and help 
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to achieve densities high enough to trigger the gravitational col- 
lapse. 

4.2. Supersonic turbulence and comparison 

The simulations with a mean Mach number of two and four 
present the same three phase structures: the ionized region, the 
shocked region, and the unperturbed region. However, in these 
cases the initial gas has more turbulent ram pressure to resist the 
ionized-gas expansion. The pressure of the ionized gas and the 
ram pressure caused by turbulence can be estimated from 



H(x -y) = 1 , if X > J or 0, if X < J. 



(3) 



pii = InuhTu 

Pturb 



(2) 



where nn and Tn ^ 7700 K are the density and the temperature 
in the ionized gas and where po and vq the density and the mean 
velocity in the initial turbulent medium. The density of the ion- 
ized gas is nearly constant in time in the three simulations, and 
does not depend on the initial turbulence, nn is around 10/cm^ 
and the corresponding pressure at t = 500 ky is given in Table [T] 
with the initial turbulent ram pressure. The ionized-gas pressure 
dominates the turbulent ram pressure for a mean Mach number 
of one and two. The consequent progression of the ionized-gas 
expansion is clearly visible in the column density plots in Fig. 
[3l However, the plots also show that this progression is nearly 
stopped for a mean Mach number of four. This is caused by the 
fact that the turbulent ram pressure dominates the ionized-gas 
pressure (see Table [T]). Another interesting difference between 
the Mach-four case and the others is the presence of many glob- 
ules. It seems that the presence of these globules in our simula- 
tions can be linked with the fact that turbulent motions dominate 
the ionized-gas expansion, this will be investigated in more de- 
tail in Sect. El 

Locally the ionization fraction changes from one in the ion- 
Table 1. Initial turbulent ram pressure and ionized-gas pressure 
at t = 500 ky 



Mean Mach number 


1 


2 


4 


Prurb/kb ait = Oky[K/cm'] 
pn/kb at t = 500 ky [K/cm^] 


3.2x10^ 
1.7x10^ 


1.2x10^ 
1.9x10^ 


4.5x10^ 
1.8x10^ 



ized region to zero in the unionized region in a very small zone 
(less than 0.01 pc) when turbulence can be neglected. Therefore 
the position at which the fraction is less than 0.5 gives the posi- 
tion of the ionization front. When turbulence is high, the position 
of the front is difficult to find since there is mixing between ion- 
ized gas and cold gas (see Fig. |4|. However global parameters 
like the mean position and apparent width of the ionizing front 
can still be studied thanks to the whole ionization fraction field. 
We used an average of the vertical profile of the ionization frac- 
tion. It transits between one in the completely ionized region and 
zero in the cold region, in between it slowly decreases because 
of the apparent width caused by the projection. We define the 
mean position of the ionizing front Zx where the average profile 
is at 0.5 and the width AZx as the distance between the 0.95 and 
0.05 positions. This corresponds to the following equations 



Z: 
AZ: 



n^-box 

x = H({X(x,y,z)h,y-0.5)dz, 
Jo 

:x = H({X),,y - 0.05) - H({X),,y - 0.95)Jz, 

Jo 



where X(x,y,z) is the ionization fraction field and H the 
Heavy-side function. Zx and AZx are shown in Fig.fTJ The front 
progression is slowed down and its width increases as the Mach 
number increases. This is caused by the bigger density contrast 
with increasing turbulence. In the three cases, the width of the 
shell is mainly caused by a projection eff'ect of the column 
density. The real shell is thin (see Fig. |4]), but is very spatially 
disperse by the inner structures of the cloud. Therefore, it will 
look wider in column density. 
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Fig. 5. 2D plot of the mass fraction at a given density and pres- 
sure (top) and at a given density and Mach number (bottom), be- 
fore ionization is turned on in the simulation with a mean Mach 
number of one. Blue to black contours are increasing mass frac- 
tion contours (blue: lE-6, green: lE-5, yellow: 5E-5, orange: lE- 
4, red: 5E-4, black: 2E-3). The dashed-red line at 4.8x10^ /cm^ 
is the maximum density achieved in a plan-parallel ID simula- 
tion for a homogeneous density at 500 /cm^. 



4.3. Ionization and thermal equilibrium 

The scheme we used (Sect. O is able to treat the ionization 
physics out of equilibrium. Therefore an interesting quantity 
to follow is the fraction of gas which is out of equilibrium. At 
equilibrium, the ionization rate compensates the recombination 
rate. From Eq. [T] the ionization fraction at equilibrium is 
X = ( -h 4/y - l)y/2 with y = o-yFy/nnlS. In the molecular 
cloud Fy is therefore X is at equilibrium. In the HII region, 
J » 1 so the ionization fraction tends to 1. At the ionization 
front, the ionization fraction can be at equilibrium and in 
between and 1, however this region is very small and can 
be neglected in the simulation (less than 1 % of the box). The 
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Fig. 4. Density cuts through the three different simulations at t = 500 ky. Left: Mach 1, middle Mach 2 and right Mach 4. 
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Fig. 6. 2D plot of the mass fraction at a given density and pres- 
sure (top) and at a given density and Mach number (bottom), 
1 My after ionization is turned on. Blue to black contours are 
increasing mass fraction contours (blue: lE-6, green: lE-5, yel- 
low: 5E-5, orange: lE-4, red: 5E-4, black: 2E-3). The dashed-red 
line is the maximum density achieved in a plan-parallel ID sim- 
ulation for a homogeneous density at 500 /cm^: 4.8x10"^ /cm^. 



thickness of the ionization front is given by 1 /nno-y which is 
of the order of 10"^ pc for ^//=500/cm^. Therefore we define 
the fraction of gas out of equilibrium as the percentage of 
cells that have a ionization fraction between 0.05 and 0.95, 
the values at diff'erent times for the three simulations are given 
in Table [J] The more turbulent the medium, the more out of 
equilibrium. It is interesting to note that this is similar to what 
was observed with the thermal in stability in the interstellar 
medium (see ISanchez- Salcedo et a l. 2002; iGazol et al.l 120051: 
lAudit & Hennebelldl2005l) . Almost one fifth of the simulation 
box at t = 1 My is out of equilibrium for a mean Mach number of 



four. This means that the hypothesis of an ionization equilibrium 
is quite bad to study the impact of ionization on a turbulent 
cloud. This is mainly caused by the ionized gas getting into a 
shadowed region. In these regions the recombination time is 
long, between 1 and 100 ky for densities around I.-IO. /cm^ and 
temperature around 8000 K as we have already explained in 
sect. [3 

The temperature is also an important quantity to monitor in 



Table 2. Percentage of cells out of ionization equilibrium. 



Time 


Mach 1 


Mach 2 


Mach 4 


250 ky 


1 % 


7% 


8% 


500 ky 


3 % 


9% 


12 % 


750 ky 


5 % 


10% 


13 % 


IMy 


6% 


12 % 


17 % 



our scheme. We solve the thermal balance between ionization 
and recombination and this allows the gas to be at states out of 
equilibrium as shown in sect. [2l This behavior can be identified 
in Fig. [8l the gas ahead of the dense parts facing the ionization 
flux is ionized but at a temperature between 2000 and 6000 K. 
This is caused by the 3D expansion of the dense ionized gas at 
the top of the structures facing the ionization radiation. Besides, 
the regions that are in the shadow of the dense structures tend to 
be overheated at almost 10 000 Kelvin whereas the equilibrium 
temperature is around 7750 Kelvin. This is caused by the fact 
that the recombination process in the shadow is done at constant 
pressure, the pressure being imposed by the hot ionized and 
exposed gas surrounding the shadowed gas. Therefore, as we 
have demonstrated in sect. [2] for the simple situation of a clump 
shadowing parts of the hot ionized gas, the temperature behind 
the clump will increase. In the turbulent simulation, the cooling 
processes for the unionized gas is also playing a role (see 
paper I for details), they will cool the recombined gas so that 
the temperature does not reach 14 000 Kelvin but equilibrates 
around 10 000 Kelvin. This result is counter-intuitive since 
one would expect that the gas stays at 7750 Kelvin during the 
recombination step, the kinetic energy of the electrons being 
radiatively lost for the system. However, the 3D versus ID 
study done in sect. [2] clearly show the importance of the lateral 
motions from the exposed gas to the shadowed gas to impose a 
recombination at constant pressure behind the structures. 
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200 400 600 800 1000 1200 1400 1600 1800 
Time [ky] 

Fig. 7. Positions of the 0.95-, 0.5- and 0.05 -transition of the xy- 
average ionization fraction. Between and the lower bound the 
ionization fraction of the gas is between 1 and 0.95, between the 
lower bound and the middle bound, between 0.95 and 0.5 and 
between the middle bound and the upper bound the ionization 
fraction is between 0.5 and 0.05. 
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Fig. 8. Temperature cut at t ^ 710 ky after ionization is turned 
on. The 3D temperature effects studied in the clump simulation 
(sect. [2] are clearly visible through the snapshot. 



5. Structures and observational signatures 

5.1. Pillars 

It is not clear how to distinguish structures on the column 
density maps in Fig. [3] because of the confusion on the line 
of sight. Therefore we will investigate some density cuts to 
identify pillars and globules. Figure [9] presents the formation 
of a pillar in the Mach 1 simulation between 200 and 700 ky. 
The pillar formed is one parsec long and presents a complex 
structure, with successive dense parts along its length. Figure 
[TOl is a map along the pillar with the direction of the velocity 
field. The structure of the field is similar to the one already 
obtained in paper I for non-turbulent set-ups. The pillar has 
two dense heads at a vertical position of 1.4 and 1.6 parsec, a 
dense base between 1.8 and 2 parsec and two holes on the side 
at 2.4 parsec. The velocity field shows that the pillar is forming 
at its base by the convergence of the dense parts of the shell 
while the holes are formed by the divergence of the dense parts 



created on a concave zone. This is the same configuration that 
was previously seen in paper I, without any turbulence. This 
suggests that the mechanism at the origin of the formation of 
pillars are not to much dependent on the turbulence. Here the 
turbulence gives only the initial conditions in terms of density 
contrasts and initial structures on which pillars will form. The 
holes on the side of the pillars in the non-turbulent set-ups of 
paper I were in fact a complete annulus around the pillars caused 
by the symmetry of our initial conditions. The major diff'erence 
with a turbulent scenario is that there is no more symmetry in 
the initial conditions. This leads to holes that are local and do 
not extend in a full annulus. 

The curvature of the shell can also be seen on spectra 



4.4 
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3.6 
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[parsec] 
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Fig. 9. Density cut at t ^ 240 ky and t ^ 710 ky after ionization 
is turned on. The white box indicates the areas of the column 
density and line-of- sight velocity spectra made in Fig. [TT] and 
[T2I The arrow indicates the direction of the line of sight; for the 
cross, the line of sight is perpendicular to the density cut plane. 

of the line-of- sight velocity. Figures [TT] and [12] show on the 
left a column density maps of the pillar and on the right the 
mass histograms as a function of the line-of-sight velocity in 
squares of 0.05 parsec. Figure [TT] is a snapshot at t 240 ky 
when the pillar is not yet formed and Fig. [121 at t 710 ky with 
the whole pillar. The line-of-sight velocity spectra are wide 
(between -2 and -\-2 km/s) before the pillar is formed. This was 
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also taking place in our turbulent simulation. 

Based on this analysis, we can infer that molecular tracers 
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Fig. 10. Cut of density field and velocity field orientation alon^ 
the pillar at t ^ 710ky. 
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Fig. 11. Left: column density on a 0.6x0.5x0.5 pc^ box around 
the area in which the pillar will form at t^^ 240 ky. Right: mass- 
weighted histogram of the line of sight velocity in the same box 
(similar to optically-thin observational line spectra). Each spec- 
trum is made on a square of 0.1x0.1 pc^ drawn on the column 
density map. The spectra are drawn between -4 and 4 km/s in 80 
bins (horizontal axis) and the mass between 10"^ and 1 solar mass 
(vertical axis in log scale). The lateral shocks can be identified 
in the wide spectra leading to a very broad averaged line width. 



already identified in paper I as the signature of the dense shell 
which is curved on itself with the two components, blue- shifted 
and red- shifted, that are going to collide to form the pillar. 
When the collision has occurred, the wide spectra is no longer 
visible and the line-of- sight velocity spectra are peaked around 
a null velocity (see Fig. El)- This demonstrates that the scenario 
identified in paper I with the study of non-turbulent set-ups is 










-J. 


-11^ 


1 1''* 


1 1 "i ' 


.J, ■ 










1 l'''l 










■± 




1 1 










Li 








1 l^^l 




JU 










1 1 1 












1 ,1 ■ 










M, 










,1., 




-X 




J, < ' 






J. 






J. 








J, 




Mill 






■ I- 




,1 




















-,Im- 






















I,- 


■J 




-,l 





Fig. 12. Left: column density on a 0.7x1.4x0.5 pc^ box around 
the pillar structure at t ^ 7 10 ky. Right: mass- weighted histogram 
of the line of sight velocity in the same box (similar to optically- 
thin observational line spectra). Each spectrum is made on a 
square of 0.1x0.1 pc^ drawn on the column density map. The 
spectra are drawn between -4 and 4 km/s in 80 bins (horizontal 
axis) and the mass between 10"^ and 1 solar mass (vertical axis 
in log scale). The lateral shocks have collided to form the pillar, 
the spectra are now peaked around km/s. 

that are optically thin should reflect the same spectral struc- 
ture as the line-of- sight velocity spectra of the simulations. 
Therefore, wide spectra could be the indication of nascent 
pillars whereas evolved pillars should present a spectra peaked 
around the velocity of the expanding bubble projected on the 
line of sight. This velocity is zero in our case since the line of 
sight is perpendicular to the direction of the expansion. 



5.2. Dense clumps at the edge of the ionized gas 

The difl'erence between the formation mechanism of pillars and 
clumps inside the shell is also clearly visible in Fig. [9l On the 
first snapshot the position on which the pillar is forming can be 
identified as the position at which the curvature of the shell is 
higher. On the other parts of the shell, the curvature is not suf- 
ficient to form pillars but will rather trigger the formation of 
clumps inside the shell. This phenomenon was already seen in 
Paper I, the higher the curvature, the longer the pillar. If the shell 
is curved enough to collapse on itself it will form a pillar. At 
high curvature, the tip of the curved shell will be stopped by the 
collapse of the matter ahead of it, which leads to the formation 
of a pillar (Right part of Fig. [131 at early times the velocity field 
of the shocked gas shows that the shell collapses on the centre 
of the initial structure). The shell around the "hill" will collapse 
quickly on the hill, with velocities that are nearly perpendicu- 
lar to the direction of the expansion. That is why the motion 
of the gas at the top of the hill is small. When the curvature is 
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low, the shell collapse on the hill but with a velocity that has 
a component in the direction of the expansion. Therefore after 
the collapse, the shocked gas keeps a velocity in the direction 
of the expansion and propagates with the rest of the shell. The 
matter is accelerated with the shell and the density increases be- 
cause the gas surrounding it converges smoothly at the centre 
(Left part of Fig. [13]). This situation was studied and compared 
in the context of RCW36 in Vela C molecular cloud (see Minier 
et al. 2012 in prep). The new observations from the Herschel 
space telescope show the presence of dense clumps around the 
HII region RCW36. Dedicated numerical simulations show that 
it is rather unlikely that these clumps were pre-existing because 
they would have triggered the formation of pillars. These dense 
clumps are located at the edge of the HII region and therefore 
they are formed by the lateral flows in the shell caused by curva- 
ture perturbations that are not high enough to form pillars. It is 
exactly the same process that is at work in the present turbulent 
simulations. From paper I, it can be estimated that the transition 
from clump to pillar formation is at a curvature radius between 
2.5 x10"^ and 5x10" ^ pc in the present situation. 

[Thompson et al.l (1201 ll) showed that the dense clumps 
around HII regions are preferentially located at the interface be- 
tween the ionized gas and the cold gas. We propose here a for- 
mation mechanism of these dense clumps in the shell thanks to 
curvature perturbations. However if the turbulent ram pressure 
dominates the ionized gas pressure, such dense clumps will not 
be present exactly at the interface, the cold gas as sufficient ki- 
netic energy to disrupt the interface and a mixing region between 
the dense shocked gas and the ionized gas is present (see Fig. 
131). Therefore the clear correlation in the observations between 
the position of the dense clumps and the interface suggests that 
in most of these regions the ionized-gas pressure dominates the 
turbulent ram pressure of the surrounding gas. 
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Fig. 13. Illustration of the eff'ect of the shell curvature on the 
formation of clumps and pillars. Density cut and velocity field 
(for the unionized gas) are taken from simulations from paper I. 
Initial curvature radius: left: 5x10"^ pc and right: 1x10"^ pc. 
Top: 20 ky after the beginning of the simulation, the shell is 
formed on the curved surface. Bottom: snapshots after 550 ky. A 
1.5-pc long pillar is formed at high curvature whereas a clump 
of the size of the initial structure is formed at low curvature. 



5.3. Globules 

Contrary to the dense clumps in the shell, globules are bubbles of 
gas disconnected from the molecular cloud and surrounding by 
ionized gas. The major diff'erence between the turbulent simula- 
tions presented in this paper and the non-turbulent set-ups stud- 
ied in paper I is the appearance of these globules in the highly 
turbulent case. As we have already discussed above, globules 
emerge because the ram pressure of turbulence dominates the 
pressure of the ionized gas. These bubbles of cold dense gas 
have enough kinetic energy to penetrate inside the HII region. 
An interesting consequence is that the motion of the globules 
is imposed by the turbulence. They could have by chance a mo- 
tion aligned with the direction of expansion of the HII region but 
most of the time, they should have a random motion direction set 
by the turbulence, and it is indeed the case in our simulation. The 
direction of expansion of the HII region is from top to bottom in 
Fig.[3l Therefore projected on the line of sight which is perpen- 
dicular, the velocity in the expansion direction is zero. The mo- 
tion of a pillar is typically aligned with the direction of expansion 
of the HII region, as it can be seen on Fig. [121 The spectra are 
centered at zero. However it is not the case for a globule. Figure 
[T4l shows a typical globule in the simulation at Mach four. The 
line-of-sight velocity spectra are systematically shifted at a ve- 
locity of -1-4 km/s which is of the order of ten times the sound 
speed of the cold gas. The shift is also highlighted in Fig. \T5\ 
The globule does not have a motion aligned with the direction of 
expansion. In a region where globules and pillars are present, the 
velocity of a pillar is the velocity of expansion of the front pro- 
jected on the line of sight, whereas the velocity of a globule will 



be the signature of the initial turbulence. The dense clumps that 
are observed in the shell also move with the expansion of the HII 
region. Therefore the same velocity shift between these clumps 
and the globules nearby can be expected. In that picture, pillars 
and clumps at the edge of the HII region are the result of a struc- 
ture dominated by the ionization dynamics whereas globules are 
turbulent dominated. This is quite diff'erent from the radiation 
driven implosion scenario, in which there is no reason for the 
globules to have a velocity diff'erent from pillars. 

This signature of the turbulent nature of a globule is sup- 
ported by recent observations in the Cygnus X r egion. Thanks to 
Hersch el Open Time and SOFIA observations (ISchneider et al.l 
I2012bh . 138 yum and CO maps are available and contain sev- 
eral pillars and globules. A pillar and a globule in this region 
were identified to be sufficiently close (and the HII region suf- 
ficiently big) to consider they should have the same direction of 
motion if the motion is only imposed by the HII region expan- 
sion. However the [CII] spectra show a line-of-sight velocity dif- 
ference of two kilometers per second, which is clearly support- 
ing a turbulent scenario for the formation of the globule against 
the radiation implosion of an isolated clump. 

The diff'erence of 2 km/s is therefore an indication of the 
level of initial turbulence. For a sonic velocity of 0.4 km/s, this 
shift indicates typically a motion at Mach 5. In our simulation, 
the globule is at Mach 10 whereas the mean Mach number of 
the simulation is 4. Therefore the velocity of the globule is an 
indication of the turbulence level but not a precise measurement. 
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Fig. 16. Probability density function of the gas for the three simulations. Left: t= ky, middle: t = 500 ky, right: t = 1 My. The first 
peak is the ionized gas and therefore is not seen in the observations. At t = 500 ky, the Mach-1 simulation presents a double peak, 
the first one corresponds to the unperturbed turbulent gas, the second one to the shocked region. In the Mach-2 simulations, the peak 
is barely visible because the turbulence is higher and because the unperturbed region still present in the box becomes small. In the 
Mach-4 simulation, the turbulence is too high, and the second peak is hidden in the peak of the turbulent cold region. 




Fig. 14. Left: column density on a 0.3x0.3x0.25 pc^ box around 
a globule in the Mach 4 simulation. Right: mass-weighted his- 
togram of the line of sight velocity in the same box (similar to 
optically-thin observational line spectra). Each spectrum is made 
on a square of 0.05x0.05 pc^ drawn on the column density map. 
The spectra are drawn between -6 and 6 km/s in 80 bins (hori- 
zontal axis) and the mass between 10"^ and 1 solar mass (vertical 
axis in log scale). All the spectra are red- shifted at a velocity of 
-h4 km/s. 
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Fig. 15. Mass fraction of the gas in the pillar and in the globule 
in function of the velocity. The bulk line-of-sight velocity of the 
pillar is at km/s whereas the bulk velocity of the globule is at 
-h4 km/s. 



5.4. Probability density function 

In this section, we investigate the density structure of the 
gas using probabil ity den s ity functions (PDFs). Numerical 
models (e.R. .Kritsuk et al.l l2007l : lAudit & Hennebelld l2010l: 
iBallesteros-Paredes et al.ll201ll) and observati ons of atomic and 
molecular gas (e.g. iKainulainen et all l2009l) show that PDFs 
are log-normal at low densities and can have more complex 
shapes at higher densities. We here plot the mass fraction 
against the density (Fig. [16]) from out simulations in a tempo- 
ral evolution and obtain that the PDF shape depends on tur- 
bulence and time. At low turbulence (Mach 1) the probabil- 
ity density function is double-peaked because of the ioniza- 
tion, similarly to the bimodal PDFs caused by the thermal in- 
stability (see ISanchez-Salcedo et al.l 120021: iGazol et al.l 120051: 
lAudit & Hennebelldl2005ir The dense compressed gas is form- 
ing a new peak in the distribution. At high turbulence, the proba- 
bility density function is shifted toward higher densities but does 
not present the double peak. This directly tell that a double peak 
is the signature of a region dominated by the ionization whereas 
a single peak is the signature that the turbulence is dominant. 
This signature was already observed in the Rosette nebula by 
[Schneider et al. (2012a^. The PDFs at the edge of the HII re- 
gions present the same double-peaked signature. It suggests that 
at the border of the bubble, the turbulence is low compared to 



the ionized-gas pressure. 

Based on the PDFs of the Mach 2 simulation, it is clear that 
the double peak disappears when there is no more unperturbed 
gas in the region on which the PDFs are done. Therefore a good 
way to look at the PDFs around a HII region is to do concentric 
PDFs centered on the central ionizing cluster. By increasing the 
radius of the region, we will first capture the peak of the shocked 
region, then a double peak should appear when the unperturbed 
region is included. Then the first peak should disappear at some 
point, when the shocked region becomes negligible compared 
to the unperturbed region. This point could indicates on which 
scale the central cluster has an impact on the global PDF of the 
cloud and therefore on which scale, the initial mass function is 
impacted by the ionization. If the turbulence is high, it is pos- 
sible that the shocked peak is hidden in the large peak of the 
turbulent cloud. In that case, no double peak will appear and this 
will be an indication that the turbulence in the region dominates 
over the ionized-gas pressure. 

6. Conclusions 

Structures at the interface between HII regions and molecular 
clouds can be classified mainly in three categories: pillars, glob- 
ules and dense clumps. Various scenarii have been investigated 
in the past to explain their formation: collect & collapse, radi- 
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ation driven implosion, shadowing effects and turbulence. We 
presented a new model to explain the formation of clumps and 
pillars in paper I based on the curvature of the dense shell formed 
by the collect processes. High curvature leads to pillars, low cur- 
vature to instabilities forming dense clumps and dips in the shell. 
This model was investigated using non-turbulent set-ups. In the 
present paper we show that the same mechanisms are at work 
in simulations with a turbulent cloud. Especially the same spec- 
tral observational signatures can be identified on the formation 
stages of the pillars. 

Furthermore we have shown that: 

- Because of turbulence, hot ionized gas can get into the 
shadow of cold dense gas, leading to a very long recombi- 
nation time for the gas. This gas can only be treated with a 
non-equilibrium model for ionization, up to 20 % of the box 
can get into that state. The equilibrium assumption is valid 
in situations in which the ram pressure of turbulence is not 
dominating the ionized-gas pressure, i.e. low turbulence lev- 
els. 

- Globules are formed preferentially when the turbulence in 
the cold gas dominates the ionized-gas pressure. Bubbles of 
dense gas have sufficient kinetic energy to penetrate into the 
HII regions forming the globules. A signature of this sce- 
nario is the line-of- sight velocity spectrum of the globules 
which is either blue- shifted or red- shifted compared to the 
spectrum of a pillar or a clump at the edge of the HII re- 
gion. This can be directly observed in regions where pillars, 
clumps and globules are present. 

- Probability density functions are double peaked when the 
turbulent ram pressure is low compared to the ionized-gas 
pressure. They could be used in observations as probes to de- 
termine the relative importance of the turbulence compared 
to the ionization. 

An important step now is to do the statistics of the forma- 
tion of dense structures at the edge of HII regions. It will be of 
great importance to conclude on the potential negative or pos- 
itive eff'ects of radiative feedback on the star-formation rates. 
Furthermore comparisons with observations will also tell how 
realistic the curvature scenario for the formation of dense clumps 
and pillars and the turbulent scenario for the formation of glob- 
ules are. 
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